#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include "vectormath.h"

inline unsigned int xorwow() {
  static unsigned int x = 123456789, y = 362436069, z = 521288629,
                      w = 88675123, v = 5783321, d = 6615241;
  unsigned int t;
  t = x ^ (x >> 2);
  x = y;
  y = z;
  z = w;
  w = v;
  v = (v ^ (v << 4)) ^ (t ^ (t << 1));
  return (d += 362437) + v;
}
inline float RandUniform() {
  const float k2Pow32Inv = 2.3283064e-10f;
  return xorwow() * k2Pow32Inv + k2Pow32Inv / 2.0f;
}

inline bool GenerateSequence(const vm::vec2 &c0, int *n, std::vector<vm::vec2> *seq) {
  vm::vec2 c = vm::vec2(0.0f, 0.0f), cn;
  *n = 0;
  for (std::size_t i = 0; i < seq->size(); ++i) {
    cn.x = c.x * c.x - c.y * c.y + c0.x;
    cn.y = 2.0f * c.x * c.y + c0.y;
    (*seq)[i] = cn;
    if (vm::dot(cn, cn) > 10.0f) {
      *n = i;
      return true;
    }
    c.x = cn.x;
    c.y = cn.y;
  }
  return false;
}
int main(int argc, char **argv) {
  float offsetx = 0.5f;
  float offsety = 0.0f;
  float zoom = 0.3f;
  int iterations = 100;
  int seqsize = 200;
  int size = 800;

  for (int i = 0; i < argc; ++i) {
    if (std::strcmp(argv[i], "-offsetx") == 0) {
      offsetx = std::atof(argv[++i]);
    } else if (std::strcmp(argv[i], "-offsety") == 0) {
      offsety = std::atof(argv[++i]);
    } else if (std::strcmp(argv[i], "-zoom") == 0) {
      zoom = std::atof(argv[++i]);
    } else if (std::strcmp(argv[i], "-iterations") == 0) {
      iterations = std::atoi(argv[++i]);
    } else if (std::strcmp(argv[i], "-seqsize") == 0) {
      seqsize = std::atoi(argv[++i]);
    } else if (std::strcmp(argv[i], "-size") == 0) {
      size = std::atoi(argv[++i]);
    }
  }
  std::printf("offsetx    = %f\n", offsetx);
  std::printf("offsety    = %f\n", offsety);
  std::printf("zoom       = %f\n", zoom);
  std::printf("iterations = %d\n", iterations);
  std::printf("seqsize    = %d\n", seqsize);
  std::printf("size       = %d\n", size);

  std::vector<unsigned int> image(size * size);
  std::vector<vm::vec2> seq(seqsize);
  for (int ii = 0; ii < iterations; ++ii) {
    for (int i = 0; i < 1000000; ++i) {
      vm::vec2 c = 6.0f * vm::vec2(RandUniform(), RandUniform()) - 3.0f;
      int n;
      if (GenerateSequence(c, &n, &seq)) {
        for (int j = 0; j < n; ++j) {
          int ix = zoom * size * (seq[j].x + offsetx) + size/2;
          int iy = zoom * size * (seq[j].y + offsety) + size/2;
          if (ix >= 0 && iy >= 0 && ix < size && iy < size) {
            image[iy * size + ix]++;
          }
        }
      }
    }
    std::printf("%.2f%%\n", static_cast<float>(ii+1) / iterations * 100.0f);
  }
  unsigned int minval, maxval = 0;
  for (int i = 0; i < size*size; ++i) {
    maxval = std::max(maxval, image[i]);
  }
  minval = maxval;
  for (int i = 0; i < size*size; ++i) {
    minval = std::min(minval, image[i]);
  }
  std::fprintf(stderr, "Density value range: %d to %d\n", minval, maxval);

  FILE *fp = std::fopen("buddhabrot.tga", "wb");
  const unsigned char tga_header[18] = {
    0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    size & 0x00ff, (size & 0xff00) >> 8,
    size & 0x00ff, (size & 0xff00) >> 8,
    32, 0
  };
  std::fwrite(tga_header, sizeof(tga_header), 1, fp);
  const float r = static_cast<float>(maxval) - minval;
  for (int y = 0; y < size; ++y) {
    for (int x = 0; x < size; ++x) {
      const int i = x + y * size;
      const int c = static_cast<int>(255.0f * (2.0f * (image[i] - minval) / r));
      unsigned char color[4] = { c, c, c, 255 };
      std::fwrite(color, sizeof(color), 1, fp);
    }
  }
  std::fclose(fp);
  return EXIT_SUCCESS;
}
